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ABSTRACT 

We study the merging history of dark matter haloes that end up in rich clusters, us- 
ing iV-body simulations of a scale-free universe. We compare the predictions of the 
extended Press & Schechter (P&S ) formalism (Bond et al. 1991; Bower 1991; Lacey 
& Cole 1993) with several conditional statistics of the proto-cluster matter: the mass 
distribution and relative abundance of progenitor haloes at different redshifts, the in- 
fall rate of progenitors within the proto-cluster, the formation redshift of the most 
massive cluster progenitor, and the accretion rates of other haloes onto it. The high 
quality of our simulations allows an unprecedented resolution in the mass range of the 
studied distributions. We also present the global mass function for the same cosmolog- 
ical model. We find that the P&S formalism and its extensions cannot simultaneously 
describe the global evolution of clustering and its evolution in a proto-cluster environ- 
ment. The best- fit P&S model for the global mass function is a poor fit to the statistics 
of cluster progenitors. This discrepancy is in the sense of underpredicting the number 
of high-mass progenitors at high redshift. Although the P&S formalism can provide a 
good qualitative description of the global evolution of hierarchical clustering, particu- 
lar attention is needed when applying the theory to the mass distribution of progenitor 
objects at high redshift. 

Key words: cosmology: theory - dark matter 



>^ ' 1 INTRODUCTION 



The Press-Schechter formalism (Press & Schechter 1974, 
hereafter P&S ) and its extensions (Bond et al., 1991; Bower 
1991; Lacey & Cole 1993, hereafter LC93) provide a nice 
framework of analytical predictions for the formation of dark 
matter haloes in cosmological models of hierarchical cluster- 
ing. Until now, these predictions have given us the best han- 
dle we possess to understand the formation of structure in 
the real universe, from sub-galactic scales up to the scale of 
galaxy clusters. It is therefore of great importance to know 
how reliable this model is, and how it fares when compared 
to numerical simulations of structure formation. 

Although the P&S predictions have been extensively 
tested against A''-body results (e.g. Efstathiou et al. 1988; 
Lacey & Cole 1994, hereafter LC94; Gelb & Bertschinger 
1994), numerical techniques are only now reaching a reso- 
lution sufficient to extend this comparison to a large range 
of masses. Due also to this limitation, most studies have 
been confined to the global mass function, while the pre- 
dicted properties of halo progenitors have been compared to 
numerical results only by LC94. 

A good understanding of the model and reliability for 
these conditional statistics is needed in order to reliably pre- 



dict several observables, like the abundance of progenitors 
of any class of astrophysical objects, their formation times, 
their merging rates, and many other properties. These distri- 
butions are particularly important for semi-analytical mod- 
els of galaxy formation (Kauffmann, White & Guiderdoni 
1993; Cole et al. 1994), where the extended P&S formalism 
is used to produce Monte Carlo realizations of the merg- 
ing histories of dark matter haloes, and these in turn are 
the starting ingredient for producing galaxy populations in 
different cosmological models. The evolution of the proper- 
ties of these populations changes substantially for different 
choices of parameters in the P&S model. 

In this paper we make a detailed analysis of the P&S 
predictions for the formation of rich galaxy clusters, and 
compare these to the properties measured in A'^-body simu- 
lations of a scale-free Einstein-de Sitter universe. The high 
mass and force resolution of our simulations is ideal to re- 
solve dark matter haloes over a large mass range. Our result 
confirms that the P&S theory can be tuned to give a good 
overall fit of the global mass function. However, we will prove 
that the P&S model best fitting the global mass function of 
our simulations is not a good fit to the statistics involving 
the progenitors of clusters found in the same simulations. 
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This result points to a structural failure of the theory, as 
the model parameter - namely the collapse threshold Sc, de- 
fined below in Section - is fixed by the match with the 
global mass function and cannot be changed when turning 
to the conditional statistics. 

In Section 2 we briefly present the A*'-body simulations 
used in our study, and describe the method used to define 
their population of dark matter haloes. In Section 3 we cali- 
brate the P&S model by fitting the global mass function of a 
set of cosmological simulations. We show that the standard 
value for the collapse threshold Sc = 1.69 provides a good fit 
to the simulations. In Sections 4 and 5 we consider a sample 
of clusters extracted from the simulations presented in Sec- 
tion 3: we present several conditional statistics of the cluster 
progenitors, and compare them with the predictions of the 
extended P&S formalism. We show that the P&S model, 
calibrated as in Section 3, is not a good description of the 
evolution of cluster progenitors. Finally, in Section 6 we dis- 
cuss our results and summarize the main conclusions of the 
paper. 



2 METHOD 

2.1 The Simulations 

The cluster A''-body simulations studied in this paper have 
been presented in Tormen, Bouchet & White (1997), where 
complete details may be found. In summary, our sample con- 
sists of nine dark matter haloes of rich galaxy clusters. These 
have been produced by tailor-made, high-resolution A'"-body 
simulations, using the initial condition resampling technique 
described in Tormen et al. (1997). The clusters correspond to 
the nine most massive objects found in a parent cosmological 
simulation (White 1994) of an Einstein-de Sitter universe, 
with scale- free power spectrum of fluctuations P{k) oc k", 
and a spectral index n = — 1, the appropriate value to mimic 
a standard cold dark matter spectrum on scales relevant to 
cluster formation. All simulations have a Hubble parameter 
Ho = 50 km s^^ Mpc^^, and are normalized to match the 
observed local abundance of galaxy clusters (White et al. 
1993). 

The parent simulation was evolved using a Particle- 
Particle-Particle-Mesh code (Efstathiou et al. 1985) with 
100"^ particles in a 256^ grid with periodic boundary con- 
ditions. The box size is 150 Mpc on a side. The cluster sim- 
ulations were evolved using a binary tree-code, in proper 
coordinates, on a sphere of diameter 150 Mpc, with vac- 
uum boundary conditions. The average cluster mass over 
our high-resolution sample is M„ ~ 1.1 x 10^^ Mq. Each 
cluster is resolved by ~ 20000 dark matter particles with an 
effective force resolution of ~ 25 kpc which is constant in 
proper coordinates. 

2.2 Identification of dark matter lumps 

We used the potential energy of particles to identify the 
centres of dark matter lumps (Efstathiou et al. 1988). We 
defined lumps by an overdensity criterion, and included all 
particles within a sphere of mean overdensity 5^ — 178, cen- 
tred on the particle with lowest potential energy. The value 
Sv corresponds to the virial overdensity in the model of a 



spherical top-hat collapse in an Einstein-de Sitter universe. 
We will call the corresponding radius the virial radius r„ of 
the lump. 

Following Tormen (1997) we limit our analysis to haloes 
with Wi, > 8 particles. We will generally call field particles all 
particles unclustered or in lumps with n„ < 8. Fig. |l| shows 
an example of the lumps found by this algorithm. The plot 
refers to a cluster in our sample, observed at redshift z = 2.8. 
Only particles which are progenitors of the final cluster are 
plotted. 



3 CALIBRATION OF THE P&S MODEL 

In this Section we calibrate the P&S model by fitting the 
predicted global mass function of dark matter haloes to the 
results of our cosmological simulations. We take the collapse 
threshold Sc as a free parameter and determine the value 
best-fitting the numerical data. We assume a top-hat filter- 
ing to define the mass variance a^{R). In this Section, and 
only here, we use two methods to define the haloes. A stan- 
dard percolation scheme - the friend-of-friends algorithm 
(hereafter FOF) - as defined in Davis et al. (1985), and a 
spherical overdensity criterion (herafter SO), as defined in 
the previous Section. We do this to demonstrate the robust- 
ness of our best fitting results. Having proven that, in the 
remainder of the paper we will only consider SO haloes. 

For the FOF we take a finking length b = 0.2. With 
this choice, groups are approximately identified by closed 
surfaces of local density p/pb ~ 1/fe"^ ~ 125 (factors of 2 or 
1/2 may apply, see Frenk et al. 1988; LC94). Following the 
notation of LC94, we will call 'FOF(0.2) haloes' the groups 
of particles found by this algorithm. For the SO haloes we 
take the virial overdensity Sv = 178, and will call them 
'SO(178) haloes'. As discussed in LC94, the algorithms we 
use to define the dark matter haloes are essentially equiva- 
lent to those used by LC94. 

3.1 Global Mass function 

The study of the global mass function requires a full cosmo- 
logical simulation. Taking advantage of the self-similarity 
of our model, we present results averaged over ten output 
times, taken from two different simulations: one is the parent 
of the cluster simulations presented above, the other is iden- 
tical to the first in the numerical parameters, and only differs 
in the realization of the initial conditions. An extensive com- 
parison between the predictions of the P&S formalism and 
the results of scale-free A'^-body simulations has been per- 
formed by LC94. At the end of this Section we will compare 
our results to theirs. 

For an Einstein-de Sitter universe with scale-free power 
spectrum P{k) oc k", the fraction of mass in haloes with 
mass M at redshift z, per interval din A/, is predicted by 
the P&S formalism to be 

dfM . a /2y/2 .M\"/2 

" ^^4-2(1?:) J' w 

where a = (ri -f- 3)/3. The time dependence enters 
only through the characteristic non-linear mass Mt{z) = 
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Figure 1. The figure shows the progenitor haloes of a cluster in our sample, seen at redshift z = 2.8; the cluster virial mass at redshift 
2 = is My = 7 X 10^^ M0. The circles show the lumps of matter found by our algorithm, based on local minima of gravitational 
potential. The radius of each circle is the virial radius of the lump. Lumps contain at least 8 particles. Overlapping circles are only due 
to projection effects. The most massive cluster progenitor is situated at position (0,0), and is formed by ^ 1300 particles. Only particles 
actually ending up in the final cluster are plotted. Therefore, progenitor haloes which carry only a small fraction of their mass into the 
final cluster have only few particles drawn within the corresponding circles. 



{'iiT/3)Rlpb{z), with pb{z) the background density of the 
universe at redshift z, and R*{z) the scale corresponding to 
a linear overdensity 5c of order unity. The specific value of 
5c, called collapse threshold, depends on the criterion used 
to identify dark matter haloes; a sensible starting value is 
5c = 1.69, the linear extrapolation at collapse for a spherical 
top-hat. We use a top-hat filter to calculate M, . The value 
of M. for Sc = 1.69 is M.|^=o = 6.16 x 10"Mq. 

The goal of this Section is to determine the value of 5c 



(that is, of M,) for which Eq. {}^ best fits the mass function 
of haloes from simulations. 

Fig. ^presents the global mass function dfu/dlnM for 
the scale-free model with n = — 1. The average over different 
output times was done as described in LC94, weighting the 
contribution of each curve in each bin by the number of 
haloes in that bin. To facilitate comparison with LC94, we 
only used haloes with at least 20 particles. 

The top left panel shows the logarithm of dfu/dlnM 
for FOF(0.2) groups (thick solid), and for SO(178) groups 
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Figure 2. Global mass function for the scale-free cosmological model with n = —1. Top Left: Thick curves indicate the values found in 
simulations, for different group finding algorithm: FOF(0.2) (solid) and SO(178) (short dash). Smooth curves show the values predicted 
by Eq. (jl]) for different values of Sc Top center: The SO(178) curve is shifted by 0.065 along the horizontal axis, showing that its shape is 
almost identical to the FOF(0.2) curve. Top right: Mass function for 10 different output times (broken curves), taken from two different 
JV-body runs of the same model. All outputs agree with the P&S prediction with Sc = 1.69 (smooth curve), and there is no systematic 
variations between the two runs or between different time outputs. This result is shown for FOF(0.2) groups, but identical results are 
obtained for SO(178) groups. Bottom Left: Logarithmic residuals from the top left panel, for FOF(0.2) groups. Error bars represent 
the Poissonian la of the mean around the observed value. Bottom Right: Corresponding residuals for SO(178) groups. Different line 
types in all panels correspond to different values of Sc. 1.3 (dotted), 1.55 (dot-dash), 1.63 (long dash) and 1.69 (solid). 



(thick dashed), together with the prediction of Eq. (|l|) for 
diflFerent values of collapse threshold: the canonical value 
5c — 1.69 (thin solid) and 5c — 1.3 (thin dotted). The choice 
of the latter value is explained below. 

From the figure we see that the mass function from 
simulations is globally flatter than the P&S prediction, espe- 
cially at the high mass end. This finding is however in general 
agreement with previous studies (e.g. Efstathiou et al. 1988, 
LC94). Despite the slight difference of shape, the A^-body 
results are reasonably fitted by the P&S prediction with the 
canonical choice 5c = 1.69, both for FOF(0.2) and SO(178) 
groups. In fact, the agreement between model and simula- 
tion curve in the mass range —1.4 ^ log]^Q(M/M«) ^ 1.4 
is better than ~ 40 percent for FOF(0.2) haloes and bet- 
ter than ~ 35 percent for SO(178) haloes. As the shape 
of the curves from A'^-body is slightly different from the 
predicted one, best fitting is a rather subjective exercise. 



However, with the choice 5c = 1.55 for FOF(0.2) haloes 
one can slightly improve the fit at high masses, up to 
logio(M/M.) ~ 1.7. 

The figure also shows that FOF(0.2) groups are slightly 
more massive than SO(178) groups. In fact, if we shift the 
SO(178) curve horizontally by 0.065 in logio(M/A/,), as in 
the top central panel, the two curves match almost exactly. 
Using this scaling between masses of haloes defined by the 
two algorithms, the best fitting value 5c — 1.55 for FOF(0.2) 
haloes correspond to 5c = 1.63 for SO(178) haloes respec- 
tively. The residuals between model and A''-body data, for 
the considered values of 5c, are plotted in the lower panels 
of Fig. |. 

The top right panel shows the 10 individual curves used 
to define the average mass function shown in the top left 
panel. Note that all outputs agree with the P&S prediction 
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with Sc = 1.69 and there is no systematic variations between 
the two runs or between different time outputs. 



3.2 Discussion 

Let us compare in more detail our results with those of LC94. 
These are based on scale-free simulations, one for each spec- 
tral index of fluctuations n — —2, n = —1, n = 0. LC94 
found that the value of 5^ best fitting the mass function of 
their simulations is Sc — 1.96 for the SO(178) haloes and 
Sc = 1.81 for the FOF(0.2) ones. With this choice, they 
could fit the A'^-body results to better than 30 per cent over 
a mass range that, for the n = — 1 simulation, corresponds to 
-1.5 <, logio M/M, <, 1.4. They also found that the canon- 
ical value Sc = 1.69 still provides a reasonable fit, with an 
accuracy better than a factor of two on the same mass range 
for the FOF(0.2) groups. 

Our results are generally in good agreement with LC94, 
as we also find that (a) SO(178) haloes are slightly smaller 
than FOF(0.2) haloes, and (b) Sc = 1.69 is an overall good fit 
to the data for groups in both identification schemes. Some 
minor differences are due to the fact that our results are 
only based on the n = — 1 model, as opposed to the three 
scale-free models used by LC94. 

Therefore, our robust conclusion is that the P&S model 
with Sc = 1.69 is an appropriate description of the clustering 
of matter in a global cosmological environment, for both 
SO(178) and FOF(0.2) dark matter haloes. In the remainder 
of the paper we will compare the predictions of this calibrated 
model to the conditional statistics of the progenitors of the 
most massive clusters found in one of the two simulations 
presented above. 

Another result is that the choice Sc — 1.3 gives clearly 
a bad fit to the global mass function (long dashed lines in 
Fig. 1^, as it predicts far too many high mass haloes. This 
result is important, as we will show below that, incidently, 
Sc — 1.3 would be the value appropriate to describe the 
conditional statistics related to the cluster progenitors. 

Finally, let us comment on the relative mass of 
FOF(0.2) and SO(178) haloes. As we mentioned before, 
haloes identified with a FOF criterion are approximately 
bounded by a surface of local overdensity p/pt « 125 for 
b = 0.2. On the other hand, SO(178) haloes are cut at a 
mean overdensity < p > / Pb = 178. Since we know the 
shape of the density profiles for these haloes (Navarro, Frenk 
& White 1996; Cole & Lacey 1996; Tormen et al. 1997), 
we can relate these two values. The virial mean overden- 
sity corresponds to a local overdensity p/pt « 10 — 30, 
weakly depending on the halo mass in terms of A/* and 
on the cosmological model. On the other hand, a local 
overdensity p/pt ~ 125 corresponds to a mean overdensity 
< p > I Ph = 500 (as found also by Frenk et al. 1988), and 
to a mass MpoFia.-i) ~ 0.75Mso{i78), with Mso(i7s) the 
virial mass. From this information, we should expect the 
FOF(0.2) groups to be less massive, not more massive than 
the SO(178) haloes. Instead we have shown above that, on 
average, MpoF(o.2) ~ 1.15M5o(i78), or about 50 per cent 
more massive than their expected mass. This difference is 
probably due to the tendency of FOF to connect structures 
along bridges and filaments, and is a quantitative evidence 
that FOF groups are, on average, very aspherical objects. 



4 EVOLUTION OF CLUSTER PROGENITORS 

We now turn our analysis to the cluster simulations. We re- 
mind that the clusters analyzed here have been obtained by 
resimulating at higher mass and spatial resolutions the nine 
most massive clusters found in one of the two A''-body runs 
presented in Section 3T . In this Section we study the mass 
function of the cluster progenitors at different redshifts, the 
evolution of the comoving density of progenitors, and their 
global infall rate as haloes merge together to form more mas- 
sive haloes. In the next Section we will concentrate on the 
most massive cluster progenitor: we will study its formation 
redshift and the time evolution of the corresponding infall 
rate. We will prove that the calibrated P&S model (which 
has Sc = 1.69), underpredicts the clustering of high-mass 
progenitors at high redshift. As mentioned above, in the re- 
mainder of this paper we only consider haloes defined by the 
spherical overdensity criterion SO(178). 

In the hierarchical clustering picture of structure for- 
mation, matter clusters on small scales first, and haloes of 
a given mass are formed by the assembly of pre-existing 
smaller haloes. Therefore, in this model, once an element 
of matter becomes part of a halo, it will be forever part 
of some halo from then on. Also, all the matter forming a 
halo comes from haloes of some mass. Finally, when a satel- 
lite lump is accreted onto a more massive halo, all of the 
satellite mass is incorporated in the more massive object. 
In this picture the formation of structure is a very lumpy 
process. This description is fairly close to what we observe 
in TV-body simulations of hierarchical clustering, although 
the actual collapse of structures is more complicated than 
its idealized version; in particular, satellites may be tidally 
stripped prior to merger with a more massive haloes, so that 
only a fraction of the satellite mass becomes part of the ac- 
creting system. Also, during violent merging events, energy 
is re-distributed among the mass elements, and some matter 
may temporarily exit the virial radius of the newly formed 
halo, or in principle even become unbound. 

Let us call, at any given redshift, progenitors of a cluster 
all haloes containing at least one particle that by 2: = will 
be part of the cluster, and also all field particles which will 
end up in the final cluster. We can ask: does all the mass 
of the progenitors end up in the final system? To answer 
we calculated the ratio Mpro{z) / of the total mass in 
progenitors (both in haloes and in field particles) at redshift 
z, to the final cluster mass. This is shown in the upper panel 
of Fig. ^. Initially, before the formation of any halo, this 
value is unity. Then, as matter starts to collapse into haloes, 
progenitor and non-progenitor particles are captured in the 
same haloes, and AIpro(z)/Mv increases up to ~ 1.2 at z ~ 1. 
It then remains at his level until z ~ 0.2, then it decreases 
to reach unity at the final time. 

What is the fate of the non-progenitor particles cap- 
tured in progenitor haloes? That is, where can we find, at 
2 = 0, the extra « 20% of matter in progenitors shown in 
the upper panel of Fig. ^? The answer is shown in the lower 
panel, where we plot the radius of spheres centered on the 
final cluster and containing different percentages of the to- 
tal progenitor mass, as identified at different redshifts. For 
example, we see that 99 per cent of the progenitor particles 
is situated, at z = 0, within roughly 2.5 virial radii, and that 
100 per cent of them is within 4.5 virial radii. This means 
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Figure 3. Top panel: fraction of the total mass in progenitors 
which will actually end up within the virial radius of the final 
cluster. Bottom panel: radius of the sphere centered on the 
final cluster and enclosing 90%, 95%, 99%, and 100% of the total 
mass in progenitors. 



that, during the cluster collapse, no particle from the pro- 
genitor haloes has actually escaped to infinity. Correspond- 
ingly, but implicitly required by the re-simulation technique 
used for the cluster simulations, no particle from infinity has 
become part of the final halo. By number, roughly 80% of 
the progenitor lumps carry at least 80% of their mass into 
the final cluster, while only 10% carry less than 20% of their 
mass, in fair agreement with the theoretical picture. 

A final remark is on resolution issues. The gravitational 
softening of our simulations is constant in proper coordi- 
nates. Since at high redshift haloes of a given mass are 
smaller, they are also more poorly resolved. As a soften- 
ing comparable to the halo size introduces a bias in the halo 
mass, in this section we conservatively consider only pro- 
genitor haloes with virial radius larger than 100 kpc, or four 
times the gravitational softening. This will result in a time 
dependent low mass cutoff in the statistics we are going to 
present. 



of haloes at redshift z, conditioned by the requirement that 
the matter will be in a cluster of given mass M2 at red- 
shift Z2 < z. For cosmological models with scale-free power 
spectrum of fluctuations: P{k) oc k", the equation reads 



dlnM 



{M,z\M2,Z2) 



exp 



1/2 fA£Y^^^Z-Z2) 



(l-a;")3/2 



Z2 



2(1 - x'^) 



(2) 



where x — M/M2 . In our case we take 2:2 = as the present 
time in Eq.(B), and M2 = M„ as the final cluster mass. 
The fraction of mass in haloes with mass larger than m„ at 
redshift z is easily obtained by integrating Eq.(^) over M, 
from rriv to infinity. The resulting cumulative mass function 
is 



/m(M > m„,z|M„ 



\V2\Mj 



(3) 



In Figs § and | we compare the differential and cumula- 
tive mass function of progenitors Eq.(|^ and Eq.(^) to the 
corresponding statistics from our simulations. Abscissa are 
normalized to the final mass of the cluster. 

Data from all nine clusters are weighted by the num- 
ber of particles in progenitors, and combined together. Each 
panel refers to a different time output, labeled by the average 
redshift over the sample. For the differential mass function, 
error bars were obtained from 100 bootstrap resampling of 
the data. That is, for each bin they indicate the la disper- 
sion of the distribution of 100 re-samples of 9 clusters each, 
each re-sample randomly drawn with replacement from the 
original collection of 9 clusters. The total number of progen- 
itors in each histogram varies from over 2000 at high redshift 
to a few hundred close to z = 0. Histograms are built using 
the virial mass of the progenitor haloes. 

It is clear from the two figures that the calibrated P&S 
model (Eqq. (||), (^), smooth solid curves) is not a good 
fit to the numerical data. Specifically, it underestimate the 
number of massive progenitors by up to one order of magni- 
tude at redshift z ^ 1. Now, it happens that the simulation 
results are reasonably fit by a P&S model with Sc = 1.3 
(corresponding to AI,\z=o = 1-35 x 10^"* Mq, dotted curve). 
However, the collapse threshold has been fixed by the fit to 
the global mass function and is not a free parameter any- 
more. Therefore, the dotted curve must be considered only 
as a rough quantitative measure of the discrepancy between 
the calibrated P&S model (which has Sc = 1.69) and the 
simulations. 

In the present analysis we have used the total mass 
of the progenitor haloes. However we found that the result 
does not change if instead one uses, for each progenitor, the 
actual fraction of halo mass ending in the final cluster. 



4.1 Conditional mass function 

In this section we fit the calibrated P&S model to the con- 
ditional mass function of cluster progenitors. In the P&S 
formalism, the mass fraction in progenitor haloes of a given 
mass can be easily calculated. Eq.(2.15) of LC93 gives 
the differential mass distribution dfM{M,z\M2,Z2)/d\nM 



4.2 Comoving density of progenitors 

Another way to look at the same results is to consider the co- 
moving density Ppro{z) of progenitor haloes in a given mass 
range [Ala, Mb], at different redshifts. This is easily obtained 
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Figure 4. Differential mass function progenitors at different redshifts. Each panel shows the mass function of the cluster progenitors, 
averaged over the cluster sample. Error bars are la of the distribution, calculated from 100 bootstrap resamples of the actual sample. 
The smooth solid curve represents the calibrated P&S prediction (Eq. with Sc = 1.69). The smooth dotted curve is a P&S model with 
Sc = 1.3. 



as an integral over Eq. (^) 



Ppro{Ma,Mt,z\M2,Z2) 



Pb 



In Mj, 



I In M, 
dfM 



dlnM 



{M,z\M2,Z2)^^i4) 



The corresponding number of progenitors per cluster in the 
same mass range is Npro{z) = Ppro{z)M2/ pb- This quantity 
is similar to the comoving density of haloes usually calcu- 
lated from the global mass function Eq. (|^), which is often 
compared to the observed abundance of some class of ex- 



tragalactic objects; it differs however in the interpretation. 
The comoving density Eq. (^) in fact refers to proto-cluster 
matter, and so it only applies to regions of the universe ini- 
tially overdense, and unusually populated. This conditional 
density can therefore be used to estimate the predicted co- 
moving density of objects which populate the environment 
of a proto-cluster. 

In Fig. ^ we show comoving density for four different 
mass ranges and for Z2 = 0. In each panel, the number of 
progenitors per comoving cubic Mpc is plotted versus red- 
shift: the symbols show the average value over the cluster 
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Figure 5. Cumulative mass function of progenitors at different redshifts. The thick solid curves are the simulation result. The thin 
smooth curves show the predictions of Eq.(b|), with 5c = 1.69 (solid line) and <5c = 1.3 (dotted line) respectively. 



sample. Virial masses for the progenitor SO(178) haloes are 
considered. 

The solid curve is the conditional density predicted 
by the calibrated P&S model. For comparison, the dashed 
line shows the corresponding global comoving density ob- 
tained by inserting the global mass distribution Eq.(^ into 
Eq. (^. The P&S model is for a final halo of mass Mi, = 
1.1x10^^ Mq , the average cluster mass for our sample. Again 
we see that the calibrated P&S model (solid curve) underes- 
timates the clustering of high-mass progenitor by up to two 
orders of magnitude at high redshift. Again, the discrep- 



ancy is roughly described by a change in collapse threshold 
{5c = 1.3, dotted curve). 

4.3 Global infall rates 

A measure of the rate at which the proto-cluster matter 
merges into more and more massive haloes is given by the 
time derivative, at any fixed mass, of the cumulative mass 
function Eq.dS). The expression is (Bower 1991) 



df{m.,,z\M^ 
dt 



d . , .dxdz 
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solid curve is the conditional calibrated P&S prediction Eq. (Q). 
The dashed line is the corresponding prediction for an uncon- 
strained environment, obtained by inserting Eq. (|^) in Eq. (Q). 
The dotted curve is a conditional P&S prediction with Sc = 1.3. 
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with xm the content of the curled brackets in Eq. More 
specifically, Eq. gives the global rate at which progenitors 
with mass smaller than rriv merge together and form pro- 
genitors of mass larger than m„, at redshift z. This quantity 
has been called infall rate by Bower (1991). Therefore, we 
keep this definition, although we will distinguish below be- 
tween the global infall rate of progenitors, given by Eq. (j^, 
and the infall rate of haloes onto the main cluster progeni- 



tor, which we calculate in Section ^j. It is in fact important 
to remember that Eq. (|^) does not contain any information 
on the relative mass of the merging haloes, nor does it dis- 
tinguish mergers onto the most massive cluster progenitor 
from other encounters between progenitors. 

It is useful to plot the information of Eq. (|^ as a func- 
tion of the lookback time from z = 0. This is done in Fig. ^ 
where we compare the predictions of Eq. (^) with the infall 
rate measured from the simulations, taking the numerical 
derivative of the curves shown in Fig. ^. Rate units are frac- 
tion of the final cluster mass per Gyr. 

Each panel shows the infall rate for a given mass thresh- 
old, in terms of the final cluster mass. The calibrated P&S 
prediction (solid curve) is systematically shifted to lower 
lookback times for ra^/Mv > 0.01, compared to the numeri- 
cal results. The dotted curve is a P&S model with 5c = 1.3, 
shown to quantify this discrepancy. 



5 FORMATION OF THE MAIN PROGENITOR 

In the previous Section we have considered the global evo- 
lution of all progenitors of a present-day cluster. A different 
way to estimate the speed at which the cluster grew is to 
look at the formation of the most massive cluster progen- 
itor, which we do in the present Section. We will focus on 
the redshift at which the main progenitor formed, and on 
the infall rate of other cluster progenitors onto it. 

Operationally, we define the main (or largest) progenitor 
of a cluster, at any given redshift, as the progenitor halo 
containing the largest fraction of the mass from the final 
cluster. In all cases this coincided with the most massive 
progenitor, although, with the present definition, the main 
cluster progenitor in A''-body simulations is not necessarily 
the most massive one, as progenitors carry different fraction 
of their mass into the final cluster. 



5.1 Formation redshift 

The cluster formation redshift is usually defined as the ear- 
liest redshift at which the largest cluster progenitor reached 
at least 50% of the cluster's final mass. The argument given 
by LC93 to estimate formation redshifts is the following. 
The probability for a cluster of present mass M2 to have 
had at least a progenitor of mass M at redshift z is simply 
the probability in Eq. (^ times the average number M2/M 
of progenitors with mass M. Since each object can have at 
most one progenitor with half or more of its proper mass, the 
probability that a cluster had a progenitor with mass XM2, 
with 0.5 < x < 1, at redshift z, is the integral quantity 

piM > xM2, z\M2,Z2) - dXr^M ^ ^ (6) 

'J ln(a:M2} 

with Z2 = 0, and df/dlnM given by Eq. (|[). For x = 0.5 



LC94 found a good agreement between this prediction and 
their A'^-body simulations, for different cosmological models 
with scale-free power spectrum. Here we compare our data 
with the same prediction, allowing however different values 
of X between 0.5 and 1, to consider a larger part of the cluster 
formation history. 

Fig. ^ shows the probability distribution function in red- 
shift for the mass of the main cluster progenitor. Symbols 
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Figure 7. The figure shows the global infall rates of cluster progenitors, as a function of redshift, for different threshold mass rriv/My. 



indicate the actual formation history of the largest progeni- 
tor, normahzed to the final cluster mass, and averaged over 
the cluster sample. Error bars give the la dispersion of the 
sample distribution. For each value of a; = M{z)/M2{z — 0), 
redshifts are indicated which correspond to the mean (thick 
solid line), and to the 5, 25, 50, 75 and 95 percentiles (dot- 
ted, dashed and thin solid lines) of the distribution obtained 
using Eq. (^. The upper panel compares the data to the cal- 
ibrated P&S model; and the lower panel to the P&S model 
with 5c = 1.3. The predictions are for a cluster with mass 
equal to the mean mass of the sample, — 1.1 x IO^^Mq. 
Symbols above the mean indicate early formation, and vicev- 
ersa. 

The figure confirms the picture above: the sample clus- 
ters form earlier than predicted by the calibrated P&S 
model. For seven out of nine objects, the mass of the largest 
progenitor is always above the mean predicted at each red- 
shift, while three clusters actually have a history at, or out- 
side, the 95% percentile of the predicted distribution. For 
comparison. Fig. |^ shows a P&S prediction with 5c — 1.3, 
which provides a better description of the numerical data. 



5.2 Accretion rates 

We now consider the evolution of the halo infall onto the 
most massive cluster progenitor. For clarity of terms, here- 
after we will call the corresponding rate the cluster accretion 
rate, not to confuse it with the infall rate discussed in the 
previous Section. 



Unlike the global infall rate, which can be calculated 
analytically, the P&S prediction for the accretion rate onto 
the main progenitor requires the use of Monte Carlo merging 
trees, as the extended P&S formalism does not contain, by 
itself, a rule on how to break up the final cluster in smaller 
progenitors as one goes back in time. To establish such a 
rule requires some extra assumptions, and different assump- 
tions lead to different ways to construct merging trees. How- 
ever, although these may lead to differences in the details of 
the merging history, they should give similar global results 
(Kauffman & White 1993). 

The expression for the infall rate Eq. (^) does not dis- 
tinguish between infall onto the most massive cluster pro- 
genitor and mergers between other progenitors. On the other 
hand, it is often useful to have the rate of infall on the main 
progenitor of the cluster, especially if one wishes to com- 
pare the latter to some observation. A classical example is 
the Butcher-Oemler effect (Butcher & Oemler 1978): the 
fraction of blue, star-forming galaxies in clusters increases 
rapidly with redshift, so that in clusters at 2: ~ 0.4 it is a 
factor of ^ 10 higher than in clusters observed locally. 

Two mechanisms have been advanced to give a physical 
explanation to the starburst of these galaxies. The first is 
the infall hypothesis (Dressier & Gunn 1983), where star 
formation is triggered by ram-pressure from the intracluster 
medium, during the first infall of an accreted galaxy in the 
cluster potential. A second idea (Lavery & Henry 1988) is 
that galaxy-galaxy collisions in the cluster environment are 
responsible for the starburst. 



© 1998 RAS, MNRAS 000, 000-000 



Assembly of Matter 11 



N 



:> 0.8 - 



N 0.6 



N 



0.4 
1 



> 0.8 



N 0.6 



0.4 



1 1 1 r 



5 = 1.69 

c 




Si 



o 
o 

CO 
Si 



0) 



CO 

S-H 



CD 

O 
O 
CO 



0.1 



0. 



0. 



0. 



2 0. 



n I I r 



"1 I I r 



"1 r 




5 10 

lookback time [Gyr] 



Figure 9. Accretion rate onto the main cluster progenitor as a 
function of lookback time. The figure shows the mass fraction 
accreted in haloes of mass m„ £ [4 X lO'^^, 10^'']Mq per Gyr by 
clusters of average final mass 6 X 10^** Mq. Squares correspond 
to the low redshift sample (L), circles to the high redshift sample 
(H). Error bars are la of the sample mean. 



Figure 8. Cluster formation redshift. Solid squares show the time 
growth of the largest cluster progenitor, averaged over the sample, 
and normalized to the final cluster mass. Error bars indicate ±la 
of the distribution. The other lines are obtained from Eq. and 
indicate the 5, 25, 50, 75 and 95 percentiles (dotted, dashed and 
thin solid) and the mean (thick solid) of the predicted formation 
redhift for a cluster with mass My = 1.1 X 10^^ Mq, the sample 
mean mass. Panel (a) compares the data to the predictions for 
the calibrated P&S model (<5c = 1.69). Panel (b) is for 6c = 1.3. 



Clusters formed in a hierarchical scenario have a his- 
tory which is qualitatively consistent with both mechanisms 
(Bower 1991; KaufTmann 1995). In fact, a cluster of given 
physical mass is a much younger object at high than at low 
redshift. Therefore, both the infall of haloes and the merg- 
ing of its progenitors are stronger at high redshift. In these 
models, the quantity used to measure the infall of galaxies 
onto the cluster has usually been the global infall rate pre- 
sented in Section 4.3 (Bower 1991), or the fraction of cluster 
mass in galaxy-size haloes (KaufTmann 1995). However, the 
accretion rate onto the main cluster progenitor is probably 
a more appropriate quantity to consider. This differs from 
the global infall rate especially early in the cluster forma- 
tion history (high lookback times), when most of the infall 
is between progenitors other than the most massive one. 

To highlight the difference between the accretion history 
of clusters observed at z = and at higher redshift, we 
divided our sample of nine objects in two subsets: a first 
one containing the five less massive clusters, which we will 
observe ai z — and will call low redshift (L) sample; and 
a second one with the four more massive clusters, which we 
will observe at 2: ~ 0.55 and will name high redshift (H) 
sample. As the average cluster mass in these two subsets 



at the time of observation is roughly the same, M„ ~ 6 x 
we are looking at clusters of the same physical size 
at both redshifts, although the (H) clusters are more massive 
in terms of M«. 

Fig. 1^ shows, for the two subsets, the evolution of the 
average accretion rate onto the main cluster progenitor, as 
a function of lookback time. Only accretion of galaxy-size 
haloes is considered: m„ G [4x lO'^^, 10^^]Mq. Squares show 
the accretion rate for low redshift (L) clusters, circles for 
(H) clusters. The corresponding global infall rate, for (L) 
clusters, is in-between that of the two upper panels of Fig. m 
For the reason explained above, the (L) accretion rate is a 
much flatter function of lookback time than the global infall 
rate, and is only slightly peaked at very high lookback times. 
On he contrary, accretion in (H) clusters mostly happens 
very recently, when it is up to four or five times bigger than 
the accretion rate for the (L) clusters. 

Fig. |l^ shows the merging tree predictions to be com- 
pared to Fig. ^. The plotted bands enclose ±lcr of the mean 
of the accretion rate on two samples of clusters with the same 
characteristics as the two A''-body subsets. Solid bands are 
the prediction of the calibrated P&S model; hatched bands 
are a P&S model with 5c = 1.3. As can be seen, accretion is 
stronger for the (H) sample, and weaker for the (L) sample, 
as in the A'^-body case. However, the predictions of different 
Sc do not differ appreciably. This is partly consequence of 
the large scatter in the predicted rates. However, the main 
reason for the similarity is that here we consider accretion of 
galaxy-size haloes, and for such haloes the two models pre- 
dict similar abundances (as shown in Fig. |^ and have similar 
global infall rates (from Fig. ^) . We therefore conclude that 
the accretion rate on the main cluster is not sensitive to the 
value of Sc. 
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Figure 10. The figure is the analog of Figure ^, but now plotting 
the prediction of the merger tree based on the extended P&S 
formalism. The bands enclose ztlcr (of the mean) around the mean 
rate. 



Finally, Fig ^ shows large variations in the accretion 
rates of the (H) sample as a function of lookback time, even 
within the same model. This result suggests that even high 
redshift clusters can sometimes exhibit, at the time of ob- 
servation, a weak accretion rate, and that this is consistent 
with the discrete nature of matter infall. Such "weakly ac- 
creting" clusters could be naturally associated to observa- 
tions of "non Butcher-Oemler" clusters. 



6 DISCUSSION AND CONCLUSIONS 

We have studied the merging history of rich clusters, using 
A'^-body simulations of an Einstein-de Sitter universe, with 
scale- free power spectrum of fluctuations, and spectral index 
n = — 1. Our results are summarized as follows. 

(i) The formation process of dark matter haloes in our 
simulations satisfies the assumptions of the hierarchical clus- 
tering model with reasonable accuracy. 

(ii) The global P&S mass function Eq. ^ provides a 
good fit to A^-body data. Although dfM{M)/d\nM from 
simulations is slightly flatter than predicted, the standard 
choice Sc — 1.69 produces a mass function which agrees 
with the A'-body data to better than 40 percent in the range 
-1.4 ^ logio(A//M.) <, 1.4, for both FOF(0.2) and SO(178) 
haloes. This conclusion agrees with earlier work (e.g. Efs- 
tathiou et al. 1988; LC94). A finer tuning in the value of 
the threshold density: 5c — 1.55 for FOF(0.2) haloes and 
5c = 1.63 for SO(178) haloes slightly extends the agreement 
at the high mass end of the mass function. 

(iii) In contrast with the previous point, the assembly of 
matter in a proto-cluster environment takes place earlier 
than predicted by the P&S model calibrated with 5c ~ 1.69. 
As a consequence, this choice underestimates the growth 



of structures especially at the high mass tail and at high 
redshift. For example, the abundance of galaxy groups at 
redshift ~ 5 is underestimated by almost two orders of mag- 
nitude. 

The discrepancy between theory and numerical data can 
be modeled by a simple variation of the collapse threshold: 
a value 5c — 1.3, which corresponds to a typical nonlin- 
ear mass M* more than two times larger than that of the 
calibrated P&S model, would give a much better fit to the 
A'-body data. However, 5c was fixed at 1.69 by the require- 
ment that the P&S model best-fits the global mass function, 
and is not a free parameter anymore. 

(iv) The accretion rate of galaxy-size haloes onto the 
main cluster progenitor is not sensitive to the choice of 5c, 
as it samples a mass range where the difference between the 
two considered models are not significative. 

The fact that the P&S model, calibrated to fit the global 
mas function, is unable to describe the clustering of a con- 
strained environment suggests that the P&S model and ex- 
tensions, in their present formulation, may be structurally 
inadequate to give a consistent description of the gravita- 
tional clustering of matter in a general context. Such a con- 
clusions would have important theoretical consequences if 
confirmed by future work. 

We stress that the present result has nothing to do with 
the different values of 5c found by several authors when fit- 
ting the global mass function (e.g. Efstathiou et al. 1988; 
Gelb & Bertschinger 1994; LC94). As clearly dicussed by 
LC94, there the discrepancies were mainly due to different 
choices for the algorithm used to define dark matter haloes 
and for the filter used to define the linear prediction. On 
the contrary, we find that the value of 5c that best-fits the 
global mass function is inadequate to fit the statistics of halo 
progenitors which come from the same simulations, and are 
defined using the same group-finding algorithm and linear 
filter. Such comparison has never been done before and high- 
lights a novel weakness of the P&S model. 

One possibile interpretation of the present results is that 
the value of the collapse threshold appropriate for the condi- 
tional statistics might depend on the mass of the final object 
relative to the value of M». In this interpretation, 5c would 
be regarded as a free parameter to be fixed as a function of 
mass. The value we found, 5c = 1.3, would then be suitable 
for rich clusters, i.e. objects with M ^ 10Af». Higher values 
might instead be required to fit the conditional statistics of 
smaller haloes, so that the standard threshold collapse value, 
5c — 1.69, would result from the average of a continuous set 
of different values over the entire halo population. 

A physical explanation of such a dependence could be 
the departure from spherical collapse (Bond & Myers 1996; 
Monaco 1995) in a mass-dependent manner; another could 
be some systematic difference in the properties of the ini- 
tial conditions for objects of different mass. To test the first 
possibility, we measured the infall pattern of matter for the 
clusters in our sample, using the axial ratios of the ellipsoid 
defined by satellite infall onto the main cluster progenitor 
(as in Tormen 1997). We did not find any systematic dif- 
ference between high and low mass clusters; however our 
sample is too small in number and in mass range to make 
a definitive statement on this point. The second possibility 
requires a more targeted analysis, again using a larger clus- 
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ter sample and spanning a larger mass range, and work is in 
progress in this direction. 

In the light of the present results, we conclude by stress- 
ing that the P&S model should be used with care when de- 
scribing the clustering properties of halo progenitors, as such 
description, in its standard formulation, is not accurate. 
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